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I. INTRODUCTION 



A. A GENERAL-PURPOSE FINITE ELEMENT PROGRAM 

A general-purpose finite element program should be able to solve a variety of 
problems from the number of disciplines: linear and non-linear, static and dynamic 
problem of elasticity, fluid mechanics, heat transfer, etc. and can solve problems of 
large size involving a variety of elements. 

A general program is going to be voluminous and complex. It is, however, 
desirable that: 

• its logic be easily understood; 

• one or many of its parts be easily modifiable; 

• it offers possibilities to tailor its facilities for the solution particular classes of 

problems. 

The program should have a modular structure, with the modules made as 
independent from one another as is practicable. The following modular operations are 
recognized: 

1. Problem Definition (data base): 

- node coordinates and element connectivities; 

- nodal element properties; 

- boundary conditions. 

2. Element Computations: 

- integration points and associated weights; 

- interpolation functions and their derivatives; 

- Jacobian matrix, inverse, and determinant; 

- element matrices and vectors: [k], [m], {f}, etc. 

3. Assembly Operations: 

- assemblage of master matrices and vectors, IK], [M], etc. 

4. Solution: 

- factorization of master matrices and solution of equations. 

5. Result: 

- output of nodal variables and other calculated quantities: gradients, 
reactions, etc. 
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Subroutines implementing the various operations described above are 
contained in all finite element codes. The flow of information between these 
operations is problem dependent; linear, non-linear, static, and dynamic problems all 
require logic of their own. 

b. mef Program 

A program of medium complexity, called MEF, implementing the techniques of 
general-purpose program that can solve a large variety of boundary value problems of 
mathematical physics. It is written in FORTRAN IV and can be easily adaptable to 
various computers. 

The main program controls the flow of all information through the functional 
blocks by transferring control to a subroutine called BLNNNN when the block calling 
card NNNN is encountered in the input file. The subroutine BLNNNN then performs 
preliminary functions such as logical unit identification, and reading of control 
parameters for the creation of various files and tables. The subroutine then calls 
subroutine EXNNNN. In all cases, subroutine BLNNNN provides appropriate default 
parameters which will be overridden by user values if specified. Subroutine EXNNNN 
then performs the major operations of the block by calling on the needed subroutines 
in the MEF library. The above protocol holds for all blocks except STOP, COMT, 
and 1MAG. All the functions of COMT and I MAG are performed by subroutine 
BLNNNN, and the function of block STOP is performed by the main program. 

The executable functional blocks contained in the MEF are: 

BLOCK FUNCTION 

SOLR Assemblage of distributed load. 

LINM Solution of linear problem with global matrix in core. 

LIND Solution of linear problem with global matrix out of core. 

NLIN Solution of stationary non-linear problem. 

TEMP Solution of unsteady problem (linear or non-linear). 

VALP Eigenvalues and eigenvectors. 

The various blocks designed for execution of the various computations have 
similar structures since they have to: 

• construct element and load matrices; 

• assemble global matrices and vectors; 

• factorize and solve the system of equations; 
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• output the results. 

Using the constructed elements and load matrices, the subroutine element library 
ELEMLB is called. This library contains subroutines that define the individual element 
types. The ELEM03 subroutine, a thirty-two node, three dimensional isoparametic 
element was developed and added to the element library. This new element allowed the 
solution oflinear elastic structures composed of homogeneous and isotropic materials. 

Multiple sample problems were developed to fully exercise use of this new 
element. An indepth investigation was then conducted to determine the limit of 
computational ability using the newly defined element to represent physical 
phenomenons. 
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II. DESCRIPTION OF THE COMPUTER PROGRAM 



A. REFERENCE ELEMENT 

To simplify the analytical expression for elements of complex shapes an element 
of reference is introduced. Such an element is defined in an abstract non-dimensional 
space with a very simple geometrical shape. The geometry of the reference element is 
then mapped into the geometry of the real element using geometrical transformation 
expression. 

A thirty-two node, three dimensional cubic element was introduced as the 
reference element. The element has eight coner nodes and twenty four mid-side nodes 
dividing each edge in three equal parts as shown in Figure 2.1. Using this reference 
element we created the fundamental matrices and vectors in subroutine called 
ELEM03, NI03, D03 and B03 to be used in element library subroutine ELEMLB of 
the MEF program. 



C 




Figure 2. 1 Reference Element. 
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B. SUBROUTINE ELEM03 

Multi-purpose program MEF is organized to contain a library of one, two, and 
three-dimensional elements for the solution of problems from a wide variety of 
disciplines. Problems from the mechanics of solids and fluids, heat transfer, etc. have 
been solved. For each element type nn, subroutine ELEMnn controls the computation 
of all matrices and vectors. 

In element type 03, subroutine ELEM03 is organized to create the fundamental 
matrices and vectors that must be numerically integrated using methods of Numerical 
Integration described in and the computations can be carried out by the following 
steps. 

1. Operation common to all element of the same type: 

- compute the weight w f and the coordinate of integration points; 

- construct the functions N (interpolation functions), the function N 
(geometrical interpolation functions) and their derivatives 

with respect to rj, £ at the points of integration. 

2. Operation for the computation of matrix [k] of each element: 

- initialize the matrix [kj; 

- for each point of integration % T : 

• construct the Jacobian matrix [J] from the derivatives with respect 
to r\, £ of function N and the nodal coordinates of the element; 

• construct the inverse of [J] and its determinant; 

• construct the derivatives of functions N with respect to x, y, z 
starting from the derivatives with respect of £,, r|, £ ; 

• construct the matrices [B] and [D]; 

• accumulate into [k] the values of [B] t [D][B]det(J)w r calculated for 
each integration point. 

3. Operation required to compute mass matrix [mj: 

- initialize [m] 

- for each integration point 2; : 

• compute Jacobian matrix and its determinant; 

• accumulate the values of {N} <N> det(J)w r into matrix [m]. 

4. Operation required to compute consistent load vectors {f}: 

- initialize {f}; 

- for each integration point \ T : 

• compute the Jacobian matrix and its determinant; 
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• accumulate the values of {N}f y det(J)w r into {f}. 

5. Operation required to compute the residue vecter {r}: 

- using the value of {f} from 4; 

- for each integration point % : 

• compute matrices [B], [D], [J] as in 2 above; 

• accumulate the product: {f} - [B] t [D][B]{u n } w r det(J) into {r}. 

6. Operation required to compute gradients <5u at points of integral: 

- for each integration point : 

• construct matrix [B] as in 2 above; 

• compute and print gradient {<5u} = [B]{u n }. 



The subroutine ELEM03 executes one operation at a time depending on the 
value of ICODE. Control variable ICODE specifies which element operation is 
desired, and expression of this variable as follows: 

ICODE = 1 initialization of the characteristic parameters of an element 
(number of nodes, number of degrees of liberty). 

ICODE = 2 operations required by a given referance element which are 
independent of the real geometry; construction of 
interpolation function N and their derivative with 
respect to ^ at the points of integration. 

ICODE = 3 construction of matrix [k] in array VKE. 

ICODE = 4 construction of tangent matrix needed for non-linear 
problems in array VKE. 

ICODE = 5 construction of mass matrix [m] in array VKE. 

ICODE = 6 computation of residual vector (r) in array VFE. 

ICODE = 7 computation load vector {f} in array VFE. 

ICODE = 8 computation and printing of gradients {<3u}. 



CODING. 



1. Evaluate coordinates, weights, functions N and their derivatives 

Integration formular for numerical integration is the following form. 



iPS 

1 = £ wjy^) (eqn 2.1) 

where 

- are the coordinates of integration point T in r^, C, system 
coordinate corresponding to weight Wj ; 
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- Wj are the weights corresponding to integration point number; 

- IPG are the total number of integration points. 

A choice of 2, 3 or 4 integration points by dimension can be made in 
subroutine GAUSS [Ref. I:p. 265], giving respectively 8, 27 and 64 integration points, 
the weights corresponding to the integration points and their coordinates. They are in 
the array called IPG, VCPG and VKPG respectively. 

Subroutine NI03 create the array VNI that contains the shape function N- 
and their derivatives with respect to ^ tj, C, system coordinate ) as Figure 2.2. 




Figure 2.2 Block Diagram for the Shape Function and their Derivative. 

2. Computation of stiffness matrix [k] of each element. 

The explicit equation of element stiffness matrix is following: 

[k] = J / J [B] t [D][B]det(J)di; dr] d£, (eqn 2.2) 

-i -i -i 



where 

- [B] is the linear strain matrix; 

- [B] 1 is the transpose matrix of [B]; 

- [D] is the matrix of elastic constants for an isotropic material. 

The stiffness matrix that has the block diagram as Figure 2.3 is the integration 
of product of the linear strain transpose matrix, the element property matrix and the 
linear strain matrix over the volume of the reference element. 
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Figure 2.3 Block Digram of Computation of Stiffness Matrix. 

In coding, we employ numerical integration formular having the following 
generic form: 

ipg 

[k] = j^w-[k° (^)] (eqn 2.3) 

where 

- IPG are the total number of integration points; 

- w- are the weighting coefficients corresponding to each integration 
point; 

- are the coordinate of the IPG integration points; 

- [k°] is the stiffness matrix at each integration point as shown in 
equation 2.4 . 



15 



[k°] = [BI^DHBJdetCJ) 



(eqn 2.4) 



Subroutine JACOB [Ref. l:p. 63], compute the Jacobian matrix, its inverse 
and its determinant. The Jacobian matrix as Figure 2.4 is obtained as the product of 

two matrices, one containing the derivatives of the geometrical transformation 

✓ 

functions with respect to the space of the reference element, and the other containing 
the real coordinates of the geometrical nodes of the element. 

< xyz > = <N(^)>[{x n } {y n } {z n }] (eqn 2.5) 

{x n }, {y n }, {z n } being the geometrical node coordinates. The Jacobian matrix is figure 
2.4 




The inverse of Jacobian matrix [J]'* as Figure 2.5 and its determinant are 
computed by the method discribed on [Ref. l:p. 44], 

Subroutine DNIDX [Ref. l:p. 64], computes the derivatives of the shape 
function N- with respect to the coordinate system of the real element using the product 
on Figure 2.6 . 

Subroutine B03 creates the matrix as Figure 2.7 that contains the strain 
components at all direction of each node in the element using output components of 
subroutine DNIDX and rearrange the new array called VBE. 

Subroutine D03 compute the stress-strain matrix (VDE) as Figure 2.S 
[Ref. 2:p. 110] for isotropic materials ( E = Young's Modulus, v = Pisson's Ratio ). 
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Figure 2.6 The derivative of shape function w.r.t. x, y, z. 

Subroutine BTDB construct the stiffness matrix by adding the product of the 
transpose matrix VBE, the stress-strain matrix VDE and the matrix VBE of every 
integration points. 

[k] = [BrtDHBJwjdetCJ) (eqn 2.6) 

3. Computation The Mass Matrix (mj. 

The element mass matrix has the block diagram as Figure 2.9 and the explicit 
equation as following: 

ill 

[m] = Jj Jj J [N] l [N]det(J) dij dr\ d$ (eqn 2.7) 

and numerical integration form is: 
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Figure 2.7 The Array VBE. 




Figure 2.8 The Array VDE. 

[m] = ^wjtNftNldetCJ) (eqn 2.8) 

i=l 

- [N] is the shape function matrix Figure 2.10 that was created by 
subroutine NI03; 

- [N] 1 is the transpose matrix of [N]; 

- wj are the weighting coefficients corresponding to each integration 
point; 

Obviously look at the product of matrices [N^andCN] as Figure 2.11 is a 
symmetric matrix. 3y convention, the mass matrix [m] contains upper half components 
and diagonal components of this matrix. 
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Figure 2.9 Block Diagram for the Mass Matrix [m]. 




Figure 2.10 The Shape Function Matrix [N]. 

4. Computation of consistent load vector (f). 

The explicit formular for the load vector {f} is: 

(0 - ( ( jWfQdeKJ) d? dl) d^ 



-1 ~1 ’I 



(eqn 2.9) 



vz 



has the block diagram as Figure 2.12 and numerical integration form is: 



vx 



ipg t 
(0 = S i w i [N] t {f vy }det(J) 



(eqn 2.10) 



vz 
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Figure 2.11 The Product of Matrices [N] 1 and [N], 

f vx , f , f yz are the force per unit volume in direction x, y, z. 

5. Computation the residue array. 

The element residual array has the block diagram as Figure 2.13 and defined 
by: 

{r} = {0 - [k]{u n } (eqn2.ll) 

where 

- {r} is the residue vector; 

- {f} is the element load vector; 

- [k] is the element stiffness matrix; 

- {u n } is the nodal values vector. 

6. Computation of the gradients and stresses at the points of integration. 

For each integration point: 

- compute strain as Figure 2.14. [Ref. 3:p. 31] 

or compacted form 

{c} = [B]{ Ui > (eqn 2.12) 
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Figure 2.12 Block Diagram for Consistent Load Vector [fj- 




Figure 2.13 Block Diagram for the Residue Array. 

- compute stress as Figure 2.15. [Ref. 3:p. 30] 
or compacted form 
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Figure 2.14 The Strain-Node Value Relation. 




Figure 2.15 The Stress-Strain Relation. 

W = IDKCj} (eqn 2.13) 

- compute coordinate of each integration points on the real element as 
figure 2.16. 
or compacted form 



{*j} = [N](x ni } (eqn 2.14) 

Using x n , y n , z n , where N = N(^, x\ v and ^ ^ are the gauss points, 

the coordinate of thirty two nodes on the real element that existed in the array 
VCORE . All the integration points are evaluated in the reference element. Using the 
product of the transfer matrix (the shape function matrix N) and the array VCORE, we 
can evaluate the coordinate of the integration points on the real element. 



22 





f > 












> 




r > 






x 




N , 


0 


( 

55 

• 

• 

• 

• 

O 


0 


0 




X 












J2 








n 




< 


Y 


> - 


0 


N , 


0 .... 0 


*32 


0 


< 


Y n 


► 




7 




0 


0 


.... 0 


0 


N 32 




z 

n i 






w - 




is 












^ * 





Figure 2.16 Evaluation of the Coordinate of the Integration Points on the Real Element. 

Note that subroutine ELEM03 executes one operation at a time depending on 
the value of ICODE. For preserving the memory locations, some of the operations 
have been performed using the same array name as VKE in both the stiffness matrix 
[k] and the mass matrix [m]. The same thing has been done with array VFE used for 
the residual vector {r} and the load vector {f}. 
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III. THE SOLUTION OF A SIMPLE PROBLEM AND DISCUSSION OF 

RESULT 



In this -chapter the preparation of input data for the computer program are 
described. A simple problem was developed and the investigation was conducted to 
determine the ability of this cubic element. 

A. ENTRY AND EXECUTION FUNCTIONAL BLOCKS 

MEF has specialized functional blocks for the entry, verification and 
organization of the data required to define a problem. Block COOR reads the nodal 
coordinates and the number of degrees of freedom of each node, it also provides 
automatic node generation. Block COND reads the boundary condition. Block PREL 
reads element properties if required for element type being used. Block SOLC reads 
the concentrated loads. Block ELEM reads the element connectivities; it also reads 
element group information when more than one element type is used, if elements have 
different properties. This block provides automatic element generation. 

Other function blocks of MEF for the execution of particular finite element 
computations use the data base constructed by entry blocks and augment it by their 
results. Block LINM assembles and solves a linear system of equations residing 
in-core. Block LIND is similar to the block LINM but the system of equations resides 
out-of-core in a mass storage device. And block STOP terminates execution of the 
problem. 

MEF provides various levels of output. The quantity of output desired from a 
given block is controlled by a parameter on the block calling card, described in detial 
[Ref. l:pp. 440-447], which ranges from 0 (the assumed value) to 4. The default value 
provides all the information needed to verify the input stream and obtain the desired 
answers while the value 1 thru 4 provide various level of verbosity. 

B. STRUCTURE MODELING 

The first step in applying a finite element solution to the problem is the selection 
of an appropriate mesh. In many case an extrapolation of the results will be required 
and hence more than one mesh will have to be selected. The mesh size solution does 
not follow any predetermined rules and will, in general, depend on the nature of the 
problem and the judgement of the analyst. An abitrary numbering convention is 
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adopted for the nodal points of the entire structure (in distinction with the numbering 
convention for the element nodes shown in Figure 3.1. A second numbering scheme is 
also needed for the elements of a mesh. 

In order to show the function of the program, as well as some observations 
affecting the use of it, the solution of a simple problem is presented in this section. 
The problem selected is that usually presented in classical texts of strength of materials 
as a cantilever beam of uniform cross section subjected to a concentrated load at the 
end of the beam. The model used for this problem is shown in Figure 3.1. Nodes 1, 2, 
3, 4, 5, 6, 7, 8, 9, 10, 11 and 12 (here arbitrarily numbered) are constrained to zero 
displacement in all directions. The concentrated load at the end of the cantilever beam 
is considered as the consistent load and are shown in Figure 3.2. 

Four different meshes which each mesh has the thickness 9, 0.9, 0.09, 0.009 
inches and has the concentrated loads 1200, 1.2, 1.2e-3 and 1.2e-6 pounds were 
employed in order to show the variation of results with mesh size. Numerical results 
for the maximum displacement at the tip of the beam are shown in Tables I and II. 




Figure 3.1 The Model used for the Problem. 
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Figure 3.2 The consistent load at the end of the beam. 

C. DISCUSSION OF RESULTS 

Elementary beam theory gives a displacement of 0.0631 inches, however this 
beam theory neglects shear deflection and three dimensional elasticity does not do so. 
Observation of the results presented here reveal that mesh refinement lead to improved 
results which eventually will converge to a certain value. The thirty two nodal points 
brick is subject to numerical bad conditionary when used with an adverse slenderness 
ratio. This ratio being defined as the ratio of the maximum element dimension over 
the minimum dimension. For example in the numerical examples treated the 
slenderness ratio for the one and eight elements representation of the four beams 
analysis are shown in the Tables I and II. Results are acceptable for the first three 
beams in each case. However a computation of approximate values of zero, first 
invariant of the element stiffness as well as condition number are produced to 
investigate the results. 

As the slenderness ratio increases the numerical of the conditioning of the 
stiffness matrices become so bad that no significant digits can be expected out of the 
solution for displacements. The same conclusion is arrived with a mesh of eight 
elements. In such a case the slenderness ratio varies from 1.6 to 1666.6. For the 
thickness of 0.09 the slenderness ratio was adequate. To reduce the ratio, more 
element must be used and tests with 20 and 40 elements were conducted. The twenty 
elements mesh gives of 0.6 to 666.6. In the use of a mesh with fourty elements the 
largest dimension become 6.0 and must therefore be reduced to 3.0 as 3.0 is the value 
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of 120/40. Runs with such a mesh were then performed to confirm our results. On the 
basis of all these runs and the value of the condition number of these matrices it seems 
that a slenderness ratio of approximately 150 could still be used to give results with 6 
significant digits in the answers, however further examples must be treated before a 
conclusion could be reached. 

These are observations essentially made on the results of a simple problem. Thus 
no firm rules regarding the use of the newly element can be established and further 
experimention with different problems has to be conducted for this purpose. 
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TABLE I 

THE DISPLACEMENT FOR THE ONE ELEMENT 



Thickness 


Slenderness 
Ratio ' 


Zero 




Zh 


^max 


^min 


Cgnril.ti.Qu. 

Number 


Displacement 


9 


13.3 


0.24e-5 


-0. 54e-5 


7. 27.11 


. 3. 72elO 


3. 14ef4 


1. 18e06 


-0.0581 


0.9 


133.3 


0. 76e-5 


-0.16e-4 


7.11.12 


1. 64ell 


7.42eH 


2. 22e09 


-0.0515 


0.09 


1333.3 


0.74e-4 


0. 16e-2 


7.11.13 


1.64«12 


7. 42.-2 


2. 2 2 el 3 


-0.0486 


0.009 


13333.3 


0. 74e-3 


0.56e-6 


7.11el4 


1.64el3 


7. 42.-3 


4. 59el5 


* 



TABLE II 

THE DISPLACEMENT FOR THE EIGHT ELEMENTS 



Thickness 


Slenderness 

Ratio 


Zero 


2>u 


Ih 


^max 


^min 


Condition 

Number 


Displacement 


9 


1.66 


0. 33e-6 


0. 23e-5 


3.21elO 


4.66e09 


1. 86e06 


2.54.03 


-0.0625 


0.9 


16.6 


0.94e-6 


0. l2e-4 


9. 12el0 


2.06el0 


1.60e04 


1.29.06 


-0.0614 


0.09 


166.6 


0.93e-5 


0.97e-4 


8. 89ell 


2. 06ell 


1. 65e01 


1. 29el0 


-0.0597 


0.009 


1666.6 


0. 93e-4 


0. 17e-2 


8. 88el2 


2. 06el2 


1.56e-2 


1.32.14 


★ 
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APPENDIX 

ELEM03 SUBPROGRAM LISTING 

The listing of subprogram ELEM03 is provided below. It includes four 
subroutines: NI03, D03, B03 and BTDB respectively. 



SUBROUTINE ELEM03( VCORE.VPRNE, VPREE. VDLE.VKE. VFE) 

*******************************x*************************************** 

32 NODES HEXAHEDRON ELEMENT FOR 3 DIMENSIONAL ELASTICITY 
EVALUATE ELEMENT INFORMATIONS ACCORDING TO ICODE VALUE 
IC0DE=1 ELEMENT PARAMETERS 

IC0DE=2 INTERPOLATION FUNCTIONS AND GAUSS COEFFICIENTS 

IC0DE=3 STIFFNESS MATRIX 

IC0DE=4 

IC0DE=5 MASS MATRIX 
IC0DE=6 RESIDUALS 
IC0DE=7 SECOND NUMBER 
IC0DE=8 EVALUATE AND PRINT STRESSES 
ELEMENT PROPERTIES 

VPREE(l) YOUNG'S MODULUS 
VPREEt 2) POISSON'S COEFFICIENT 
VPREE(3) SPECIFIC MASS 

********************************************************************** 

IMPLICIT REAL*8(A-H ,0-Z) 

COMMON /COOR/NDIM 
COMMON /ASSE/NSYM 

COMMON /RGDT/IEL, ITPE. ITPE1 , IGRE, IDLE.ICE, IPRNE, IPREE, INEL, IDEG, 

1 ipg icod£, IDLE d , ineLo , ip6o 

COMMON /ES/M a MR,MP 

DIMENSION VC6RE( 1) ,VPRNE( 1) ,VPREE( 1) , VDLE( 1) , VKE( 1) , VFE( 1) 



CHARACTERISTIC DIMENSIONS OF THE ELEMENT 




DIMENSION VNIX(96) ,VNI(3 



GKED(3) 



— DIMENSION OF MATRIX D, NUMBER OF G. P. 

DATA IMATD/6/ . IPGKED/3 ,3 ,3/ 

DATA ZERO/O. 060/,DEUX/^. 60/,X05/0. 5D0/,RADN/. 57295779513082302/ 
DATA EPS/1. D-6/ 

DATA NNNNNN/0/ 



CHOOSE FUNCTION TO BE EXECUTED 



GOTO (100, 200, 300, 400, 500, 600, 700, 800), ICODE 
100 IDLE0=96 
INEL0=32 
IPG0=27 
RETURN 

*********************************************************************** 

EVALUATE COORDINATES .WEIGHTS , FUNCTIONS N AND THEIR 
DERIVATIVES AT G.P. 

*********************************************************************** 

200 CALL GAUSS( IPGKED.NDIM, VKPG, VCPG , IPG) 

IF(M. LT.2) GOTO 2^0 
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WRITEfMP ,2000) IPG 

2000 FORMAT ( / 1 5 , ' GAUSS POINTS ' /10X , ' VCPG ' ,25X, 1 VKPG ' ) 

10=1 

DO 210 IG=1 . IPG 
I1=I0+NDIM-1 

WRITE(MP .2010) VCPG( IG) , ( VKPG( I ) , 1=10,11) 

210 I0=I0+NDiM 

2010 FORMAT! IX, FI 3. 9, 5X. 3 FI 3. 9) 

220 CALL ,N 103 ( VKPG , VNI ) 

IF(M. LT.2) RETGRN 
I1-4*INEL*IPG 

WRITEIMP ,2020) (VNI(I) , 1=1 , II) 

2020 FORMA I (/iX .FUNCTION AlsID DERIVATIVES'/! IX ,8E12. 5)) 

RETURN " 

** ******************************************************************** 

EVALUATE ELEMENT STIFFNESS MATRIX 

********************************************************************** 

INITIALIZE VKE 

300 N I N 1=4656 

IF(NSYM. NE. 0) NINI=9216 
DO 310 1=1 .NlNI 
310 VKE(I)=ZERG 

--- FORM MATRIX D 

CALL D03! VPREE ,VDE) 

LOOP OVER THE G. P. 

1 1= 1+ 1 NE L 
DO 330 IG=1 , IPG 

EVALUATE tHE JACOBIAN, ITS INVERSE AND ITS DETERMINANT 

CALL JACOB fVNinn .VCOftE,NDIM,INEL,VJ,VJl,DETJ) 

C=VCPG( IG)*DETJ 
DO 320 1=1.36 
320 VDE1( I )=VDEfl )*C 

PERFORM MATRIX B 

CALL DNIDX ( VNIfl 1) , VJ1 ,NDIM, INEL, VNIX) 

CALL B03( VNIX , INEL, VBE) 

CALL BTDB ( VKO, VBE, VDE 1 .IDLE, IMATD.NSYM) 

330 I1=I1+4*INEL 
RETURN 

400 CONTINUE 
RETURN 

********************************************************************** 
EVALUATE THE MASS MATRIX 

********************************************************************** 

00 NINI=4656 

IF(NSYM. NE. 0) NINI=9216 
DO 510 1=1 ,NINI 
510 VKE(I)=ZER0 

LOOP OVER THE G. P. 

I D IM1=NDIM— 1 
I DEC L=C ND I M+l ) * I NE L 
I 1=1+INEL 
12=0 

DO 550 IG=1 , IPG 

CALL JACOB( VNI( 1 1) , VCORE.NDIM, INEL, VJ.VJl.DETJ) 

D=VCPG( IG)*DETJ*VPREE( LL) 

ACCUMULATE MASS TERMS 

IDL=0 

DO 540 J=1 , INEL 
JJ=I2+J 

J0=l+IDL*(IDL+l)/2 
DO 530 1=1 , J 
11=12+1 

C=VNI( II)*VNI(JJ)*D 



30 



ooo ooo o o ooo ooo ooooo 



VKE(J0)=VKE(J0)+C 
Jl=J0+IDL+2 
DO 520 11=1, IDIM1 
VKE(J1)=VKE(J1)+C 
520 Jl=Jl+IDL+3 
530 J0=J0+NDIM 
540 IDL=IDL+NDIM 
I 1=1 1+IDECL 
550 I2-I2+IDECL 
RETURN 

******************************************************** 

EVALUATE THE ELEMENT RESIDUAL 
******************************************************** 

— form MATRIX D 

600 CALL D03(VPREE, VDE) 

INITIALIZE THE RESIDUAL VECTOR 

DO 610 ID=1 , IDLE 
610 VFE(ID)=ZERO 

LOOP OVER THE G. P. 

I 1=1+1 NEL 
DO 640 IG=1 . IPG 

EVALUATE THE j A C0BIAN 

CALL JAC0B(vNI(I1) a VC0RE,NDIM.INEL.VJ,VJ1,DETJ) 

EVALUATE FUNCTIONS D(Nl)/D(X) 

CALL DN IDX( VN 1(11) , VJ1 ,NDIM, INEL, VNIX) 

EVALUATE STRAINS AND STRESSES 

EPSX=ZER0 

EPSY=ZER0 

EPSZ=ZER0 

GAMXY=ZERO 

GAMYZ=ZERO 

GAMZX=ZERO 

ID=1 

DO 620 IN=1 , INEL 
UN=VDLE( ID) 

VN=VDLE(ID+1) 

WN=VDLE( ID+2) 

C1=VNIX[ IN) 

IN1=IN+INEL 
C2=VNIX( INI 
IN2=IN1+INE 
C3=VNIX(IN2 
EPSX=EPSX+C1*UN 
EPSY=EPSY+C2*VN 
EPSZ=EPSZ+C3*WN 
GAMXY=GAMXY+C1*VN+C2*UN 
GAMYZ=GAMYZ+C2*WN+C3*VN 
GAMZX=GAMZX+C1*WN+C3*UN 
620 ID=ID+3 

C1=VCPG( IG)*DETJ 
C2=VDE(2 )*C1 
C3=VDE(22)*C1 
Cl=VDE(irCl 

SIGX=Cl*EPSX+C2*EPSY+C2*EPSZ 
SIGY=C2*EPSX+C1*EPSY+C2*EPSZ 
SIGZ=C2*EPSX+C2*EPSY+C1*EPSZ 
TAUXY=C3*GAMXY 
TAUYZ=C3*GAMYZ 
TAUZX=C3*GAMZX 

FORM THE RESIDUAL 



* *** * 
* *** * 



* * * * * 
* * * * * 
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C3=VNIX( IN2) 

VFE( ID)=VFEfID)+Cl*SIGX+C2*TAUXY+C3*TAUZX 
VFE( ID+1 )=VFE(ID+1)+C2*SIGY+C1*TAUXY+C3*TAUYZ 
VFE(ID+2)=VFE(ID+2)+C3*SIGZ+C2*TAUYZ+C1*TAUZX 
630 ID=ID+3 
640 I1=I1+4*INEL 
RETURN 

********************************************************************** 

EVALUATE BODY FORCES , FX , FY.FZ PER UNIT VOLUME 
(FOR GRAVITY FX=0 ’ FY=0 , RZ=-VPREE(3) 

********************************************************************** 

700 FX=ZER0 
FY=ZER0 
LL=3 

FZ=-VPREE(LL) 

DO 710 1=1.96 
710 VFE( I)=ZER6 

IDECL=(NDIM+1)*INEL 
DO 730 IG=1, IPG 

CALL JAC0B( ^NICIl+INEL) , VCORE.NDIM , INEL, VJ , VJ1 , DETJ) 
DX=VCPG(IG)*DETJ 
DY=DX*FY 
DZ=DX*FZ 
DX=DX*FX 
12=11 
13=1 

DO 720 IN=1, INEL 

VFE( I3)=VFE( I3)+DX*VNI( 12) 

VFE( I3+1)=VFE( I3+1)+DY*VNI( 12) 

VFE( 1 3+2 )=VFE( 1 3+2 )+DZ*VN 1(12) 

12 = 12+1 
720 13=13+3 

730 1 1=1 1+ 1 DEC L 

RETURN 

******************************************************************* 
EVALUATE AND PRINT STRESS AT G. P. 

******************************************************************* 
800 WRITEfMP , 2080) IEL 

2080 F0RMAT(//' STRESSES IN ELEMENT' ,15/' P.G. ,6X, 'X' , 1 IX , Y ,11X, 

1 ' Z ' ,9X, 1 SIGX 1 ,8X, SIGY' ,8X, SIG2 ,7X, ‘ TAUXY ' 7X TAUYZ' ,7X, 
2'TAOZX'/) 

FORM THE MATRIX D 

CALL D03(VPREE, VDE) 

LOOP OVER THE G. P. 

I1=1+INEL 

12=0 

DO 820 IG=1 , IPG 

EVALUATE THE JACOBIAN 

CALL JACOBf VN I f 1 1} . VCORE , NDIM. INEL , VJ , VJ 1 , DETJ ) 

EVALUATE FUNCTIONS D(NI)/6(X) 

CALL DNIDX(VNI( II) , VJ1 ,NDIM, INEL, VNIX) 

COMPUTE STRAINS AND COORDINATE AT G. P. 

EPSX=ZER0 

EPSY=ZERO 

EPSZ=ZER0 

GAMXY=ZERO 
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GAMYZ=ZERO 

GAMZX=ZERO 

X=ZER0 

Y=ZER0 

Z=ZER0 

ID=1 

DO 810 IN-l.INEL 
UN=VDLE(ID) 

VN=VDLE( ID+1) 

WN=VDLE(.ID+2) 

XN=VCORE(ID) 

YN=VC0RE(ID+1) 

ZN=VC0RE(ID+2) 

C1=VNIX( IN) 

IN1=IN+INEL 
C2=VNIX( INI) 
IN2=IN1+INEL 
C3=VNIX( IN2) 

IN1=IN+I2 
C4=VNI( INI) 
EPSX=EPSX+tl*UN 
EPSY=EPSY+C2*VN 
EPSZ=EPSZ+C3*WN 
GAMXY=GAMXY+C1*VN+C2*UN 
GAMYZ=GAMYZ+C2*WN+C3*VN 
GAMZX=GAMZX+C1*WN+C3*UN 
X=X+C4*XN 
Y=Y+C4*YN 
Z=Z+C4*ZN 
810 ID=ID+3 



COMPUTE THE STRESSES 



SIGX=VDE(1)*EPSX+VDE(2)*EPSY+VDE(2)*EPSZ 

SIGY=VDE(2)*EPSX+VDE(1)*EPSY+VDE(2)*EPSZ 

SIGZ=VDE(.2l*EPSX+VDE(2)*EPSY+VDE( 1)*EPSZ 

TAUXY=VDE(22)*GAMXY 

TAUYZ=VDE(22)*GAMYZ 

TAUZX=VDE( 22 j*GAMZX 

WRITEfMP ,2090) IG,X,Y,Z,SIGX,SIGY,SIGZ,TAUXY,TAUYZ,TAUZX 
2090 F0RMAT( 1# ,15, 9E12. 5 ) 

I2=I2+4*lflEL 
820 1 1=1 1+4*INEL 

RETURN 
END 

SUBROUTINE NI03(VKPG, VNI) 

*********************************************************************** 

TO EVALUATE THE SHAPE FUNCTIONS N AND THEIR DERIVATIVES W. R. T. 
KSI , ETA, DZETA 
INPUT 

VKPG(NN)= COORDINATES OF POINTS IN KSI , ETA, DZETA 
OUTPUT 

VNI = THE SHAPE FUNCTION N AND 

THE DERIVATIVE OF SHAPE FUNCTION W. R. T. KSI , ETA, DZETA 
*********************************************************************** 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION VN 1(3456) , dORRF(32,3) .VKPG(81) ,VN(32,3) ,RN(32) 

DATA ((CORRF(I.J), 0=1.3). 1=1 32) 

* /-l. do'-l.Dd.-l. DO,-. 333333333333333,-1. D0.-1. DO, 0. 3333 

*33333333333,-1. D6 -1. Dd.l. D0,-1. DO -1. D0,1. DO.-. 3333^33333^3333 . 
*-l. DO , 1. DO,. 33333^33333^333 -1. D0.1. DO.l. D0,-i. D0 A . 333333333333^33 
*1. D0.-1. D(l.-. 3333333333333^3,1. Dd.-l. do -1. DO. 1. dO.-l. DO.-l. DO, 
*d. 333^33333533333,-1. D0,-1. DO -. 335333335333335,-1. dO,-l. dO A -l. dO A 
*-. 333333333333333 ! 1. DO.-l. D0 A -. 333333333333333,1. D0.1. DO.-. 5333335 
*33333333,-1. D0.1. do- 533333533333333.-1. DO.-l. D0 A . 533335333333333 
*,1. DO.-l. DO a . 353333533333333.1. D0 A 1. Dd A . 333533333533333 ,-1. DO.l. DO 
**. 333533333533333,-1. DO.-l. Dd, 1. Dd.-. 353333333333333 .-1. DO .1. dO A 
*. 333333333333333 a -1. D0 a 1. D0,1. DO , -1. DO.l. DO.l. DO,-. 353333353333533 
* , 1. DO , 1. DO , . 333353333353333 ’ 1. DO 1. DO, I. DO , 1. DO ,. 333333333333333 , 
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c 

Cl 

C99 



*1. DO. 1. DO.-. 333333333333333.1. DO. 1. DO -1. D0 A 1. DO, 1. D0,-1. DO,. 3333 
*33335333353, l.DO.-l. DO,-. 335333353333533,1.60/ 



DO 1 1=1.32 

WRITE?*, $9) I ,(C0RRF( I ,J) ,J=1 ,3) 
FORMAT ?2X,I2 3(2X,E12.6)) 



JJ=1 
DO 10 NN= 



,81,3 



X=VKPG(NN) 

Y=VKPG(NN+1) 

Z=VKPG(NN+2) 

PRINT*'X.Y,Z 
At Conner 

NODE # : 1,4,7,10,21,24,27,30 

DO 100 1=1.2 
II=20*(I-l)+l 
IT=I 1+9 

DO 100 J=I I « IT , 3 
Xl=C0RRF(J,r 
Y1=C0RRF(J 2 
Zl=C0RRFfJ 3^ 

CF1=9. D0/64. 60 



.. T ... DO/9. DO 
RNfJ)=CFl*(l. D0+X*X1)*( 1 
*“■' D0+Y*Yl)* 



100 



. J,1)=CF1 

* +2 

VN(J ,2)=CF1* 
* + 2 . 

^VN(J ,3)=CF1* 



)o*x: 



. D0+Y*Y1)*(1 
(1.D0+Z*Z1)* 



D0+Z*Z1)*(-CF2+X*X+Y*Y+Z* 
(Xl*(-CF2+3. D0*X*X+Y*Y+Z*Z 



V 



1. D0+X*X1)*(1. D0+Z*Zl)*(Yl*(-CF2+X*X+3. D0*Y*Y+Z*Z) 
)0*Y) 

1. D0+X*X1)*(1.D0+Y*Y1)*(Z1*(-CF2+X*X+Y*Y+3.D0*Z*Z) 



AT MIDSIDE 

NODE § : 2,3,8,9,22,23,28,29 



DO 200 K=1 .2 
IK=20*(K-l) 

DO 200 1=1,2 
II=6*(I-l)+2+IK 
IT=I 1+1 

DO 200 J=I I , IT 
X1=C0RRF( J , 1 
Y1=C0RRF( J , 2 
Z1=C0RRF( J ’ 3, 
CF1=81. D0/S4. DO 
CF2=1. DO/9. DO 
RN(J)=CF1*C 



200 



VN 

VN 

VN 



1. D0-X*X)*(CF2+X*X1)*(1,D0+Y*Y1} 
*(1. D0+Y*Y1)*( 1. D0+Z*Zl)*(Xl-2. D 
*Y1*£1. D0-X*X)*(CF2+X*Xlr(l.D0+ 



‘J , l)=CFl 

;j,2)=CFl _ .. _ - 

J 3)=CF1*Z1*(1. D0~X*X)*(CF2+X*X1 



(1. D0+Z*Z1) 
0*X*CF2-3. D0*X*X*X1) 
D0+Z*Z1 j 



)*f l! D0+Y*Y1 



AT MIDSIDE 

NODE § : 5,6,11,12,25,26,31,32 



DO 300 K=l,2 
IK=20*(K-lJ 
DO 300 1=1,2 
II=6*(I-l)+5+IK 
IT=I 1+1 
DO 300 J=II. 

X1=C0RRF(J , 1 
Y1=C0RRF( J , 2 . 

Zl=C0RRFc J , 3/ 
RN(J)=CFWi.D 0+X*X1)*(1 
J , 1 )=CFl*Xl*(l. D0-Y*Y 



300 



VN 

VN 

VN 



:j 2 = 

0 , 3 * 



=CF1*( 1. D0+X*X1) 
=CF1*Z1*(1. D0+X^ 



D0-Y* 

*( CF2 
1. D0+Z*Z1 
)*( 1. D0-Y 



Y)*(CF2+Y*Y1)*(1. [ 
+Y*h)*(i.D0+z*zi: 
' 1-2. D0*Y*( 



Y)*(CF2+Y*Y1) 



0+Z*Zl) 

F2-3. D0*Y*Y*Y1) 



AT MIDSIDE 

NODE # : 13,14,15,16,17,18,19,20 
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DO 400 J=13,20 
X1=C0RRF( J , 1) 
Y1=C0RRF(J ,2) 
Z1=C0RRF(J,3^ 



RN(J)=CFl*{ 1. D0+X*X1)*( 1. D0+Y*Y1)*(1. D0-Z*Z)*(CF2+Z*Z1) 
VN(J,1)=CF1*X1*( 1. D0+Y*Y1)*( 1. D0-Z*Z)*(CF2+Z*Z1) 



. ,_,1)=CF1*X1*(1.D0+Y*Y1)*(1. D0- 
VN(J 2)=CFl*Yl*f 1. D0+X*XlWl. D0- 
400 VN(j’3)=CFl*fl. D0+X*X1)*(1. D0+Y*Y 
DO 4 I 0 L= 1.3^ 



Z*ZWCF2+Z*Zi; 

Yl)*(Zl-2. D0*Z*CF2-3. D0*Z*Z*Z1) 



VNI(JJ)= 



|=RN( I 



410 JJ=JJ+ 

DO 420 K=1 ,3 
DO 430 1=1 32 
VNI(JJ)=VN( I , K) 

430 JJ=JJ+1 
420 CONTINUE 
10 CONTINUE 
RETURN 
END 

SUBROUTINE D03( VPREE.VDE) 

*****************************&***************************************** 

TO FORM MATRIX D ( 3 DIMESIONAL ELASTICITY) 

INPUT 

VPREE = ELEMENT PROPERTY 

VPREE(l) = YOUNG'S MODULUS 
VPREE(2) = POISSON'S RATIO 

OUTPUT 

VDE = MATRIX D 

*********************************************************************** 

IMPLICIT REAL*8 (A-H.O-Z: 

DIMENSION VPREECi; , , 

v. ..... .. QEUX/2. 0D0/ 



DATA ZERO/0. 0D07 
E=VPREE( 1) 

LL=2 
X= ~ 

Cl 
C2 



VPREE(LL) 
=E*(UN-X)/(( 
=C1*X/(UN-X 
UN-r 



(A-H.O-Z) 

,iWi E £>do},i 

UN+X)*(UN-DEUX*X)) 



C3=C1*(UN-DEUX*X)/(DEUX*(UN-X)) 

DO 10 J=l,36 
10 VDE(J)=ZEftO 
VDE(1)=C1 
VDE(2)=C2 
VDE(3)=C2 
VDE( 7 )=C2 
VDE(8)=C1 
VDE(9)=C2 
VDE( 1 3 )=C2 
VDE( 14)=C2 
VDE(15)=C1 
VDE( 22 j=C3 
VDE(29)=C3 
VDE(36)=C3 
RETURN 
END 

SUBROUTINE B03(VNIX , INEL, VBE) 

*********************************** 

TO FORM MATRIX B (3 DIMEMSIONAL ELASTICITY) 

INPUT 

VNIX = DERIVATIVES OF SHAPE FUNCTION W. R. T. X,Y,Z 
OUTPUT 

VBE = MATRIX B 

*********************************************************************** 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION VNIX(32,1),VBE(6,1) 

DO 10 1=1,6 
DO 10 J=1 ’96 
10 VBE( I , J)=6. 0D0 



************************************** 
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c 

c 



FORMATION OF MATRIX B 



DO 20 1=1,32 
11=3*1-2 
12 = 11+1 
13=12+1 
KK=1 



C1=VNIX(I ,KK) 
C2=VNIXfI (KK+l)) 
- (KK+2)) 



********************** 



C3= - , X . A -- 

VBE( 1 , Il)=Cl 
VBE(2 I2)=C2 
VBE(3 I3)=C3 
VBE(4 Ili=C2 
VBE(4 I2)=C1 
VBE(5 I2)=C3 
VBE(5 I3)=C2 
VBEf 6 I1)=C3 
20 VBEf 6 I3J=C1 
RETURN 
END 

SUBROUTINE BTDBfVKE, VBE. VDE , IDLE , IMATD , NSYM) 
************************************************* 

TO ADD THE PRODUCT B(T).D.B TO THE ELEMENT MATRIX K 
INPUT 

VBE = MATRIX B 
VDE = MATRIX D 

OUTPUT 

VKE = ELEMENT MATRIX K 

*********************************************************************** 
IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION VKE( 1),VBE( IMATD ,1),VDE( IMATD ,1),T(576) 

DATA ZERO/O. 0D07 
IJ=1 

IMAX=IDLE 
DO 40 J=1 , IDLE 
DO 20 11=1 , IMATD 
C=ZERO 

DO 10 J 1=1, IMATD 
C=C+VDE( II J1)*VBE(J1 ,J) 

T(I1)=C 

I F( NSYM. EQ. 0) IMAX=J 
DO 40 1=1 , IMAX 
C=ZERO 

DO 30 Jl=l, IMATD 
C=C+VBE(J1 ’ I)*T(J1) 
yKE^IJJ=VK^(IJ)+C 

RETURN 
END 
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